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A POSTERIORI ERROR ESTIMATES FOR DISCONTINUOUS GALERKIN 
METHODS USING NON-POLYNOMIAL BASIS FUNCTIONS. 

PART I: SECOND ORDER LINEAR PDF 

Lin Lin^ and Benjamin Stamm ^ 


Abstract. We present the first systematic work for deriving a posteriori error estimates for 
general non-polynomial basis functions in an interior penalty discontinuous Galerkin (DG) 
formulation for solving second order linear PDEs. Our residual type upper and lower bound 
error estimates measure the error in the energy norm. The main merit of our method is 
that the method is parameter-free, in the sense that all but one solution-dependent constants 
appearing in the upper and lower bound estimates are explicitly computable by solving local 
eigenvalue problems, and the only non-computable constant can be reasonably approximated 
by a computable one without affecting the overall effectiveness of the estimates in practice. As 
a side product of our formulation, the penalty parameter in the interior penalty formulation 
can be automatically determined as well. We develop an efficient numerical procedure to 
compute the error estimators. Numerical results for a variety of problems in ID and 2D 
demonstrate that both the upper bound and lower bound are effective. 

1991 Mathematics Subject Classification. 65J10, 65N15, 65N30. 


1. Introduction 

Let O be a bounded domain. We consider the development of a posteriori error estimates for the 
following second order linear PDE 


— Alt -h Vu = /, 


in O, 


( 1 ) 


using the discontinuous Galerkin (DG) formulation with general non-polynomial basis sets. 

Such equation arises in many scientific and engineering problems such as in electromagnetism, geo¬ 
physics, quantum physics, to name a few. In order to solve Eq. Q in practice, it is desirable to reduce 
the number of degrees of freedom for discretizing Eq. Q to have a smaller algebraic problem to solve. 
While standard polynomial basis functions can approach a complete basis set and is versatile enough 
to represent almost any function of interest, the resulting number of degrees of freedom is usually large 
even when high order polynomials are used. Non-polynomial basis functions are therefore often em¬ 
ployed to reduce the number of degrees of freedom, and are widely used to solve Eq. Q and other 


equations, including the planewave basis set for solving Helmholtz equation 12,26 , the heterogeneous 
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multiscale method (HMM) and the multiscale finite element method 13 for solving multiscale ellip¬ 
tic equations, and the various non-polynomial basis set used in quantum chemistry such as the Gaussian 
basis set 


10 , atomic orbital basis set 16 , and adaptive local basis set 20 etc. 


Besides solving the equation, it is also often desirable to assess the accuracy of the numerical solution 
via a posteriori error estimates and to design approximation spaces that result in a uniform distribution 
of the error in space to achieve best accuracy for a given number of degrees of freedom. In this paper 
we focus on the a posteriori error estimates in the interior penalty DG formulation 21 27 . 


The DG formulation has the advantage that it formally relaxes the continuity constraint of basis 
functions at the inter-element boundary, and is therefore particularly suitable for incorporating general 
basis functions, which are difficult to match at the inter-element boundary. 


1.1. Previous work 

Compared to the many existing works on a posteriori error estimates using polynomial basis functions 
in the DG formulation [T^[l7 23 , it is much more difficult to develop systematic a posteriori error esti¬ 


mates for general non-polynomial basis functions. One of the important reasons is that approximation 
and scaling properties of the function space spanned by non-polynomial basis functions, which are key 
to a posteriori error estimates, are generally difficult to deduce. For instance, Amara et al developed 
the upper bound error estimates for the Helmholtz equation in planewave basis enriched DG method, 
and the error is measured in the L^-norm. Kaye et al 18 developed the upper bound error estimates 


for solving linear eigenvalue problems using non-polynomial basis functions in a DG framework, which 
generalizes the work of Giani et al 11 for polynomial basis functions. However, the assumption of 


approximation properties on the function space is in general difficult to verify. Though not in the DG 
framework, Ohlberger 22 developed the a posteriori error estimates for the HMM method for elliptic 


homogenization problems. 

The difficulties of a posteriori error analysis for general non-polynomial basis functions are largely 
due to the lack of credible methods for measuring the ratio of the error using different norms, defined 
in proper function spaces. For instance, approximately speaking, in a residual based error estimator, 
the constants associated with the residual requires the estimation of ratio of the error measured using 
L^-norm and the i7^-norm. The scaling properties of such constants with respect to the increase of 
the number of basis functions on a particular element can be rather intrigue for non-polynomial basis 
functions. The estimation of such constants is already complicated for polynomial basis functions or 
planewave basis functions, not to mention the case when the non-polynomial basis functions come from 
numerical solution without a analytic recipe, or even worse, the basis functions do not in practice form 
a complete basis set with only saturating accuracy. 


1.2. Contribution 

To the extent of our knowledge, this is the first systematic work for deriving a posteriori error 
estimates for general non-polynomial basis functions in a DG framework. Our upper and lower bound 
error estimates are residual type estimators for the error in the energy norm. In our formulation, all 
but one basis-dependent constants appearing in the upper and lower bound estimates are explicitly 
computable by solving local eigenvalue problems. For solution with sufficient regularity (for instance 
u G the only non-computable constant can be reasonably approximated by a computable one 

without affecting the overall effectiveness of the estimates. While the requirement of regularity 

appears to be a formal drawback in the context of a posteriori error estimates, the main goal of this 
work is to develop a posteriori error estimates for general basis sets rather than for /i-refinement, and the 
difficulty of general basis sets holds even if the solution has regularity. Therefore we think our 

method can have important practical values. As a side product of our results, the penalty parameter in 
the interior penalty formulation is also automatically computed, and the computed constants guarantees 
that the coercivity of the resulting DG bilinear form for Poisson’s equation. 

We develop an efficient numerical procedure to compute these constants. Both the formulation and 
the practical implementation of our method are independent of how the basis functions are generated. 
Although the numerical procedure is developed for general non-polynomial basis functions, we find that 
the procedure, when applied to standard polynomial basis functions, generates constants are even more 
accurate than the analytical asymptotic result. Numerical results for a variety of problems in ID and 2D 
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indicate that both the upper bound and lower bound are sharp, and the effectiveness of the estimators 
holds even at the level of each element. 

1.3. Outline 

The rest of the paper is organized as follows. After an introduction to some technical results in 
section]^ we start with the derivation of the upper bound a posteriori error estimates for the Poisson’s 
equation in section without the potential term V. We then generalize the derivation of the upper 
bound error estimates to indefinite problems with the potential term, as well as the lower bound error 
estimates in section]^ We elaborate in sectionon the numerical methods for computing the constants 
appearing in the upper and lower bound estimates needed in our analysis. Finally, we present numerical 
results in section before we conclude in section followed by an appendix. 

2. Preliminary results 


2.1. Mesh, broken spaces, jump and average operators 

Let n = (0,1)'^, d = 1, 2,3 and let /C be a regular partition of ft into elements k G 1C. That is, we 
assume that the interior of k H k', for any k, k' G /C, is either an element of K., a common face, edge, 
vertex of the partition or the empty set. For simplicity, we identify the boundary of P in a periodical 
manner. That means, that we also assume the partition to be regular across the boundary dCl. We 
remark that although the assumption of a rectangular domain with periodic boundary condition appears 
to be restrictive, such setup already directly finds its application in important areas such as quantum 
chemistry and materials science. However, the analysis below is not restricted to equations with periodic 
boundary condition. Other boundary conditions, such as Dirichlet or Neumann boundary conditions 
can be employed as well with minor modihcation. Generalization to non-rectangular domain does not 
introduce conceptual difficulties either, but may lead to changes in numerical schemes for estimating 
relevant constants in section if the tensorial structure of the grid points is not preserved. 

Let N = {Nk)i.^^jc denote the vector of the local number of degrees of freedom on each element 
K C 1C. Let Vat Y]s[{k) by any piecewise discontinuous approximation space on a partition 1C 

of the domain Cl. It is important to highlight that little is known about the a priori information of Vjv 
except that we assume that each Yn{k) contains constant functions and that Yn{k) C (k), so that 
the traces of Yvn on the boundary Ok are well-defined for all vn G Yn{k), for all k G 1C. We denote 
by H^{k) the standard Sobolev space of L^(K)-functions such that all partial derivatives of order s G N 
or less lie as well in L‘^{k). By H^{1C), we denote the set of piecewise iL'*-functions defined by 

LG(/C) = {u G I uU G H^{K)yK G 1C} , 

also referred to as the broken Sobolev space. We denote by H}^{Cl) the space of periodic iL^-functions 
on H. We further define the element-wise resp. face-wise scalar-products and norms as 

= ^(u, w)„ and ||u||k; = (u, u)^. 

kGK. 


The L^-norm on k and Cl are denoted by |j • ||k and || • ||o, respectively. 

The jump and average operators on a face F = k H k' are defined in a standard manner by 


ivj = + and H = -I- 

|Vu} = 5(Vu|«,-b and |Vu] = + VuU'n„/, 

where denotes the exterior unit normal of the element k. Further we state some standard results 
Finally we recall the standard result of piecewise integration by parts formula that will be employed 
several times in the upcoming analysis. 

Lemma 2.1. Let v,w G H'^{1C). Then, there holds 


[(Au,ui)„-b (Vu, Vu;)k 


5 yz + (Vf'. H)a« 
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2.2. Projections 

For any element k G JC, let us denote by IIq : L^{k) -g K the L^(K)-projection onto constant functions 
defined by 

(npW, = (w, Vw G K, 

that is explicitly given by IlgU = dx. On H^{k) we define the following scalar product and norm 

(r, w)*,„ = (IIqU, now)^ + (Vu, Vu;)„, (2) 

1 

IMU.k = 

for all v,w G H^{k) and the corresponding projection 11^ : H^{k) —>■ Yn{k) by 

{II%v,wn)*,k = iv,WN)i,,K VrcAT G V7 v(k). (3) 

Then, it is easy to see that this projection satisfies the following properties 


(?; - n^?;,c)K = 0, 

Vc G K, V?; G H^{k), 


or equivalently expressed as nQ(?; — n^w) = 0. 

This implies that 


(V(?; - n^?;), ViCAr)^ = 0, 

ywN G VAr(K),V?; G H^{k), 

(4) 

||V(?;-n^?;)|U< ||V?;|U, 

Vu G H\k), 

(5) 

lk-n^i^lk.« < lkll*.«> 

Vv G H\k). 



2.3. Local scaling constants 


In this section, we are going to define some local constants that will be used in the upcoming a 
posteriori error analysis. We start with defining the local trace inverse inequality constant for each 
K G JC defined by 


d 


K 


sup 

n{i^) 


|j V VjSf 
IkjvlU.K 


> 0 . 


Further, let 


_ P' U 

S-K = sup T—r- 

vGH^(k), IwIU.k 

v1.Yn(k) 


and 


sup 

v±.Yn{k) 


Mkii 

IklU,. 


where _L is in the sense of the scalar product (•, •)*,« defined by ([^. 


Remark 2.2 (The computation of the constants a^, and d,.j). We provide more details in 
on how these local constants can be approximated by solving local eigenvalue problems. 


Lemma 2.3. Let k G K., v G H^{k). Then, there hods that 


Ik - llVulU, 

||u - n^ulla^ < ||Vu||«;. 


Proof. The proof consists of simply combining the definition of resp. b^ and the stability of the 
projection described in ([^ 


II?; - n^?;|U < Ik - n^?;|U,K = a„ ||V(?; - n^?;)|U < ||V?;||k, 

since nQ(?; — n^u) = 0. The proof for the second inequality is almost identical. 


□ 
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3. Poisson’s equation 

As has been motivated in the introduction we start with a simple model problem that however reflects 
the difficulties associated to the discontinuous Galerkin method using non-polynomial functions. The 
problem then reads: find u G n such that 

- Am = /, in O, (6) 


for some / G L^( 0 ). 

Given a piecewise constant and positive penalty function 7 such that = 7k G ffi"*" for all k G JC, 
the discontinuous bilinear form is defined by 

a{w,v) = Y^ [(Vm;, Vm)k - KVm;, [mDok - + ^(M> > 


for any w,v G H^{JC) and for 6* G K. Note that this is equivalent to the somewhat more standard 
notation 


a{w,v) = {Vw,Vv)k - ({Vm;}, |M])jr - 6 »(|m;], |Vm})^ -h (tfM. H)j 


with 7 f = fy} and where denotes the face-wise L^-inner product over all faces of the mesh. 

Note that the choice of 9 = 1,-1 corresponds to the symmetric and non-symmetric interior penalty 
discontinuous Galerkin (SIPG [^27 or NIPG [^) method, respectively. The former case results in a 
symmetric bilinear form. 

Then, the discontinuous Galerkin approximation is defined by: Find mjv G Vat such that 


a{uN,VN) = if, vn)q, Vmat G Vat. 


(7) 


In this context we define the following broken energy norm by 


IMP = E [ 11 ^^ 11 ' +^IINIIL 


Vm G H\IC). 


( 8 ) 


Observe that 



kG/C 


with 


I|Vm||2 + ^||[m]|| 


2 

dK'> 


and that this is indeed a norm as 7 > 0. As usual, the penalty parameter 7 needs to be chosen carefully 
to ensure coercivity. Even when polynomial basis functions are used, the choice of an optimal 7 is not 
completely trivial and related discussions can be found in [T][^ . The scaling in the element sizes and the 
polynomial order is however known 15 25 . The involved constants are resulting from applying trace 


and inverse inequalities, but no inverse inequality is known if general non-polynomial basis functions 
are employed. To have a precise idea of the values of the combined trace and inverse inequalities for 
the generic non-polynomial basis functions spanning Vat, we propose here to use the local constants 
that were defined in Section In consequence, we can give a precise value of the piecewise constant 
function 7 that is needed to ensure coercivity of the bilinear form a(-, •). This is stated in the following 
lemma. 


Lemma 3.1. Under the assumption that ||| • ||| is a norm (which is assumed here since 7 > Oj and if 
additionally 7k > | (1 + 9)"^ (dK)^ for each n G 1C, then the bilinear form is coercive on Yn, i.e., there 
holds 


< a(MAr, Mat), Vmat G Vat. 
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Proof. Since for any G Vat we have Vfjv = '^{vn — HgLiAr) and ||i;Ar — IlgWAilU.K; = IIVtatUk we can 
develop 

a{vN,VN) = 5 X! (l + ^')(V(r’w-noi;Ar),|i;Arl)a/^+7K||bw]||L 

kGK. 

- 5 X! (l + 6»)d«||i;Ar-noi;Ar||*,^||[wjvl||a/^+7K||bjv]||L 

kGK. 

> 5 E + (^« - II 

kGK. 

and obtain 

||||i^iv||P < a{vN,VN) 

for the particular choice 7 k > 5(1 + d)^ (dK)^- Note however that for the particular choice of 0 = —1, 
7 k stills needs to be positive in order that ||| • ||| is indeed a norm. □ 

3.1. Error representation 

Define the scaled error function ip = |||^~^”||| and develop 

|||m - Ujvl = E “ UN),yp)K + ^{{u - Ujv], M)aK 

kGK 

= a{u - utv, (/?) + E 1“ “ watDok- 
kGK 


We prefer to work with the scaled error function tp for sake of a simple presentation of the upcoming 
error analysis. Observe that due to the regularity of u G which implies |m] = 0, and since u is 

indeed the solution of § , there holds 

a{u,p) = E [(Vu, V</5 )k - 5 (Vu, M)aK = -iAu,p)n = if,p)Q. 

kGK 

On the other hand, since G Vat is the DG-solution solution of 0 , we obtain 

-a{uN, p) = -a{uNj p - Pn) - (/, PN)n, 

for any p^ G Vat- Thus, using the integration by parts and Lemma [2.1[ we can develop 

-a{uN,p) = E [“ (Vuat, V(93- Pn))k + KVuat, {p - PnDok + f (|uAr], V(</5 - PN))dK 




- ^(Iww], -‘^wDok -{f,PN)n 


E [(AuAr,V5 - Pn)k - 5([VuAf],‘^ - PN)dK + f (luAll, v((^ - PN))dK 


kG/C 


- ^([mw], - (^wDok -{f,PN)Q, 


and obtain the error representation eqnation 

|||u-UAr||| = E [if + Aun,P-Pn)k - ^{lyuNl,p-pN)dK 


( 9 ) 




- 7K([MAr], ip-PN)nif)dK - ^([matI, Vy)+6»V</J7v)9K 
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3.2. A posteriori error estimation 

After recalling that we assumed that u S we start by introducing the constant d“(uAr) defined 

by 


d“(MAr) 


||V(u - UN)-ni^\\oK. 
||V(m - UAr)|U 


and define the constant c„ by 

= d“(MAr) + d«;|6»|. 

We note that in practice, the constant d“(Mjv) can not be evaluated since u is unknown. See in the 
upcoming numerical examples how we deal with this term. 

Remark 3.2. Observe that d“(itjv) is bounded by the constant 


|1V('U - VN)-n^\\o^ 

sup 

«ivGVjv(K) II V(m - UAr)||„ 


< oo, 


which, in turn, is independent of the approximation un (but still depends on the exact solution u and 
the approximation space Yn)- 

Define the following estimators 


?7r,k = a^ll/ + AuatIIk, 

(10) 

?7f,« = ^IIIVuArlllo^, 

(11) 

= (bK, “f ■^)II |mm| II Ok, 

(12) 


in order to state the first Theorem. 

Theorem 3.3. Let u G n H^{IC) be the solution of Q and un GYn the DG-approximation 

defined by 0. Then, we have the following a posteriori upper bound 

|||r - Mm III < I ^ [bR.K + + bJ.K 

\kG/C 



Proof. Using the triangle inequality, observe that 

||(V(^ + 6»V(/jAr)-n^||aK < |lV(/5-n„||a„ + |6»|||VvJM-M„||aK < d“(MAr)||V vj||«, + d«,|6 *|||V(/jjvIU- 

So far, the results where valid for any arbitrary discrete function G Y]\[. In this proof we consider 
the particular choice pn\k = so that we can easily state 

l|V(/jivlU < IIV(^IU 

by splitting p = Tlfjp + {p — and using the orthogonality relation Then, there holds 


||(Vv? + 6'V(pAr)-n«,||a„ < (d“(wAr) + d^l^l) ||V(/5||„ < Ck||V(p||k, 

^^-^ 

— Ck 


(13) 


by applying a simple triangle inequality. 
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If we appl y th e Cauchy-Schwarz inequality to the error representation formula Q, in combination 
with Lemma |2.3l Eq. (131 and another Cauchy-Schwarz inequality, we have (recall that (p = 

Ih - Mjvl < ^ [11/ -f AujvlU llv? - V^atIU + ^IIIVuAr]||a«; \w - + 7K||[wAflllaK \W - ^nWok 


k,GK. 


5llI^«Arlll9«; II(VV5 + 0'^PN)-nK\\aK 


< 


11 / + AmatIIk + ^|||VuAr]||p„ -I- + ^)|||MArll|9K II 


k,GK. 


< 


[a«, 11/ -f AmatII^ -f ^lllVnArllloK + 


Ok 


\k^K. 




+ + Vj,k, 


\k^K. 


□ 


4. Second order indefinite problems 

In this section we consider the more general indefinite equation: find u G iL^(O) such that 

— Au + Vu = f, in O, (14) 

for some / € L^(0) and where we only assume that V G L°°(fl) is bounded and that the operator 
—A + V has no zero eigenvalue. For the particular choice of E = —k^ G K this framework includes the 
Helmholtz equation. The DG-bilinear form is provided by 

a(w,v) = [(Vic, V'u)„ -f {Vw,v)^ - ^(Vtc, |z;])aK - |([■u;], V'u)^^ -f ^(|w], , 

such that the DG-approximation is defined by: Find un G Yn such that 


a{uN,VN) = {f ,VN)n, ^vnGYn, 


(15) 


and we keep the definition of the broken energy norm of (§. Of course the choice 7 „ =2(1-1- \9\Y (d^)^ 
does not imply coercivity of the bilinear form in this setting any more. We assume that 7 ^ has been 
chosen by the user to insure that DG-problem has a unique solution and focus on how to quantify the 
error a posteriori. Observe that whenever the DG-problem is not uniquely solvable, the solver of the 
numerical system typically reveals the lack of well-posedness. The following analysis requires that the 
DG-solution satisfies (15). 


4.1. Computable upper bounds 

We first introduce a modified norm. For this consider V+ and Y defined by 14 = max(l/ 0) > 0 and 
V- = max(—F, 0) > 0 so that V = V+ — V and \V\ = V+ +V-. Then, define 


Ilk 11^ = Y 

kGK. 


\Ml = \\Yv\\l + \\V^^v\\l + ^\M\\l,, 


Vn G H\]C). 


Applying similar arguments as in Section the following error representation can be developed 

Ik - unW = Y “ VUN, P-Pn)k + {V-{u- Un), p)k 

kGK. 

~ ^Y + Tfcdwjvl, + iluNl,Yip+eYpN)aK ■ 

kGK. 


(16) 
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Redefining the residual as 


= a^ll/ + Aun - VunWk, 


the following bound can be developed. 


(17) 


Theorem 4.1. Let u G H^(Q) n be the solution of (14| and um G Vjv the DG-approximation 

defined by (151. Then, we have the following a posteriori upper bound 


u — Un 



+ Vf.k 



\\V-hu-u^)\\l 
|||m - unW 


where defined by 0 and riY^K,V],K defined by ([TT|)"([T^. 

Proof. This estimate can be obtained by applying the Cauchy-Schwarz inequality to the error represen¬ 
tation equation (161 similar as in the proof of Theorem 3.3 Only the additional term 


{V-{u-un),t)k = 


\VJ{u-unWk 


\\u - ujvll 


is not estimated. □ 

||U_^ (u—■ujv)||^ 

Remark 4.2. For V- G L°°, the term — '|||„_^^||| — is small compared to the upper bound estimator 
in the limit of complete basis set. On the other hand, when only a small number of basis functions are 
used, this term can become large, and the upper bound error estimator can underestimate the true error 
in energy norm. 

4.2. Computable lower bounds 

The goal of this section is to derive computable lower bounds of the approximation error. We note 
that the following theory applies also to the Poisson problem with V = 0. 

Observe that 

^J,k; — (h^ “b ) II I^A^J II dn — \J (b^ ~2 ) III a WAI III K • 

Second, for any face F of Ok, denote by k' the adjacent element such that F = dnOdn' such that there 
holds 

= xllI^^JvlIlL = xllI^(^-'“w)lllL ^ H (ll^(«-«Af)U •»^«IIf + ||V(u-Mjv)U/ • n«/|||.y 

FGdu 

(18) 


Recall that uj{k) is the patch consisting of k and its adjacent elements sharing one face, then 


17L < ^ H K'(«a')II^(^-«A^)IU')^ ^ d“,(ujv)) ll^(“-«^)ll 

\ Ki ^ OJ V ^ ' 


k' ^uj{k) 


k' 


Further, let g^^ be a smooth non-negative bubble function with sup,^g^ ^^(a;) = 1 and local support, 
i.e. supp((7k) C k, which in turn implies that PkIok = 0. Finally, let us denote the residual by 
R = f + Aun — Vun and define 


o-K 


mu 

hi mi 


Denote by G Hhh the solution to 


= Vg^^R, 


on K, 
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SO that 


VK,K = a.K.\\R\\K = crK\\giR\\l = (TK / 5 k 

J K 

= -a^ A{u- un) g^R- 


= 


V(u — Wat) • V (5k R) — V{u 


< cr„||V(M - uiv)IUI|V(5k R - V5k)||k, 


A{u — un) + V{u 
■ Un) 

- Un) ■ V(/3 k 



R 


and in consequence 




F-'^^ivlU 


< crK||V(5K-R- (Pk) 


The results above indicate that 

|||?x — wa/-|||k > max 


^R,k, 

Or K Oj K 


iM-MArL(K) > 


5f,k 

Cf,k 


where, denoting by |w(k)| the cardinality of the set uj{k), we use the definitions 


L E iiv<. + yiiH 


■ \uj{k)\ 

' ^ k'Gw(k) 


laK) 


(19) 


and 


||i?|U||V(&Ki?-(pK)|U 

CR.K - 1 /n 5 

\\bi^ R\\l 

Cf,k = max d“,(wAr), 

k'GCjJ{k) 

Cj,K — \ ^ ' 


We summarize the results in the following proposition. 


Proposition 4.3 (Local lower bound). Let u S H^{Ll) n H^{JC) be the solution of (14) and un & Vjv 
the DG-approximation defined by (15). Then, the quantity 


t- I Vr,k 5f,k 5j,k I 

U = max<( -,-,-^ , 

^R,k ^F,k ^J,k 


is a local lower bound of the local error 

max {IIIu-M at Ik, |m - UAr|a;(K)} • 

Remark 4.4. Since in practice, the nominator as well as the denominator of any of those fractions 
might become very small, these ratios are not numerically stable. It turns out that 

^ _ 5r.k + 5f,k + 5j,k 
Cr,k + Cf,k + Cj^K 

is numerically more robust and still meaningful as it replaces the maximum by the average. 


On a global level, the following result holds. 
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Proposition 4.5 (Global lower bound). Let u € be the solution of (14) and un & Vjv 

the DG-approximation defined by (15). Then, there holds that 






GK 


VR,k. + Vf,k. + V2,k. 


where 


(cr,^ + 

b'W = = max (f + 


< |||u - Un\\ 


FedK 


FedK 


Proof. Observe that by (18) there holds 


'III 


^ XI IE(«-“Af)U-»^«lll + X ^ X IE(w-UAr)U' -n^ 

kGK kG/C F'gOk /^G/C F^Ok 

^X^ X IE(^^-«w)U-’^«lll + X X Vll^(“-“Ar)U •ri/.lll 

kG/C F^Ok kGK F^Ok 


X X (¥ + ¥) 


K • iI'kIIf ~ ^ ^ ^ ^ { 
/*G/C FgOk 


I V (zi '^iv) I K, * '^K. II F 


k^K F^Ok 

kGK. kGK 


Then, using the other local estimates for ? 7 r_k and given by (19) yields 

X ^ 3 X + 'iD ^ 3 X (cL + b^(«)d“(wAr)^ + |||u - Mjvl 

kG/C kGK. 


kGK. 


□ 


5. Practical strategies for estimating the constants 

In this section we discuss how to compute the constants dK,a,.j,bK as defined in Section in the 
a posteriori error estimator for general non-polynomial basis functions in the discontinuous Galerkin 
framework. The basic strategy is to discretize the infinite dimensional representative space H^(k) using 
a finite dimensional space such as high order polynomials, and to replace the various inner products 
defined in Section]^ by discrete bilinear forms using Gauss quadrature. With the help of these bilinear 
forms, d^, a„,b„ can be estimated by solving an eigenvalue problem, locally on each element k. 

5.1. Finite dimensional discretization 

For simplicity let k = [0, h]‘^, d= 1,2,3 and all quantities be real. We start the discussion with d = 1, 
i.e. K = [0, h]. All numerical quadrature are to be performed using the Legendre-Gauss-Lobatto (LGL) 
quadrature with Ng points. The LGL grid points are denoted by {yj}^fi, and the corresponding LGL 
weights by The Lobatto quadrature implies that 

yi = 0, VN, = h, 

which facilitates the description of the boundary integrals as in the estimate of d^ and b^. The LGL 
grid points {yj}ffi correspond to a unique set of Lagrange polynomials of degree (Ng — 1), denoted by 
{Pjix)}^li, and satisfy 


Pj{yt) = Sij, l<i,j<Ng, 
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where 6ij is the Kronecker <5 function. We can then approximate v G H^(k) using the linear combination 
of Lagrange polynomials as 

v{x) « 

i=i 

The sequence of spaces of polynomials of degree Ng being dense in H^(k) implies that, for any 
V G H^{k) and any e > 0, there exists Ng and ^ Pat^ such that 


II K 


< e 


and similarly 




< £, 


if choosing Ng large enough. That is, elements in H^{k) can be approximated, in the sense of and 
with any desired accuracy by elements of PiVg- This motivates us to work in Pat^ instead of H^{k) 
for Ng large enough. We assume that Ng is large enough so that the above approximation error in 
the local and iL^-norms are very small. Further, for functions u,v G PAig, the LGL quadrature for 
computing the inner product (m,u)k converges rapidly with respect to the increase of Ng. 

We denote by u = (ui,..., vn^)’^ the column vector corresponding to the coefficients of u G PAig, and 
denote by T = (j/i,... ,yNg)'^, w = (wi,... ,uj]s[g)'^ the column vector corresponding to the LGL grid 
points and weights, respectively. With a slight abuse of notation we can compute the inner product 
using linear algebra notation as 

Ng 

{u,v)k, = ''^^UjUjjVj = u^Wv, ( 20 ) 

1=1 

where W = diag[r(;] is a diagonal matrix with the entries of vector w on the diagonal entries. 

The Lagrange polynomials also induce a differentiation matrix D of size Ng x Ng, defined as 


Dii=p'j{Vi), l<i,j<Ng. 


( 21 ) 


Taking the derivative of a polynomial yields 


Ng 

V'{x) = '^p'g{x)Vj. 

1 = 1 

Let v' = ... ,v'{yN ))^ be the column vector of the derivative quantity v'{x) on the LGL grid 

points, then 

v' = Dv. (22) 

Eq. (221 shows that the differentiation matrix maps the values of a function to the values of its derivative 
on the LGL grid points. Using the differentiation matrix, inner products of the form {u',v')k can be 
expressed in linear algebra notation as 

(u', u')« = (DufWiDv) = u^{D^WD)v. (23) 

In order to compute the inner product (u, u)*^k we also need to compute (nQU,!!))!!)^. Note that 

noli = 7^(1,u)« = 


with IkI = h. Then 


(ngU,ngu)K = t—^u^ww'^v\k\ = ( w-Nw^ ) V. 


Therefore the inner product {u, u)*_k can be computed as 


(u, u)*_K = « ( D^WD + Wy-w^ ) V. 


(24) 





TITLE WILL BE SET BY THE PUBLISHER 


13 


We also need to compute inner products on the boundary dn. In ID, v\dK.ix) is completely described 
by two points i;(0) and v{h), which is given by the discretization on the LGL grid points as vi and 
VNg- Define the weight vector at 0-dimension as w = (1, 0,..., 0,1)^, and W = diag[ui], then the inner 
product on the boundary can be expressed as 


Similarly 


iu,v)aK = UiVi + UNgVNg =u^Wv. 


{u', v')dK, = -I- u'jq = u^D^WDv. 


(25) 

(26) 


The inner products (201, (23), (24) and (26) are sufficient for estimating dK,aK,b„ for d= 1. 


Now we generalize all the definition above to d > 1. Though in practice we only consider d = 2,3, 
the formalism developed here holds for any dimension. For any x G k = [0,d]‘^, we denote by a: = 
..., , with being the component of x along the Z-th dimension. Then the set of Ng LGL 

grid points in the dimension d is given by 


= {Viu-da = (%i> ■ • • < Ji. • ■ • .Jd < Ng}- 

We define the tensor product of d matrices ..., of size Ng x Ng as 


(27) 


i (0 


1 . . . — TT A' 

— 11^1,Ji’ 
1=1 


1 < *1, jl, ■ ■ • Gd, jd < Ng 


which can be written in a compact form as 


(28) 


a = 6^a^^\ 


(29) 


From the computational point of view, it is more convenient to rewrite the tensor product A as a matrix 
by stacking the ii,... ,id and ji,..., indices, respectively. In other words, we can view A as a large 
matrix of size Ng x and each matrix element corresponds to a matrix element Axj, 

with the index 

d d 

I = i + Y^Y- i)fvd-i), j = 1 + 

Note that when d = 2, the stacked representation of the tensor product of A^) and A^^^ is the Kronecker 
product of and A''^\ 

We also dehne a special case for the tensor product of d vectors ..., of size Ng. By viewing 
each as a matrix of size Ng x 1, we have 


1=1 


(i) 


1 < Jl, ■ , jd < Ng. 


(30) 


Eq. (30) can be written in a compact form as 


v = (^v 
1=1 


(0 


(31) 


By stacking the indices ji,... ,jd together, we can view u as a vector of size N^ and each element 

corresponds to an element vj with = 1 -I- J2'i=iUi ~ l)-^g* Using the notation of tensor 
product, the set of LGL weights is described by a vector 


1=1 


w' 


(32) 







14 


TITLE WILL BE SET BY THE PUBLISHER 


Similar to the ID case, each LGL grid point uniquely corresponds to a Lagrange polynomial 

d 

Ph,-.jdix) = 

It can be readily seen that 


;=i 


Pi 




1=1 


As in the ID case, a polynomial u{x) defined on k can be expressed using the Lagrange polynomials as 
U{X)= Ph,-Jdix)uiy3i,-dd) = Y Pjl,-Jdix)Uj^,...,33- (33) 

^<jl,- -,jd<Ng i<jl,- -,jd<Ng 

Denote by = diag[u>[‘^l] as a matrix of size Ng x Ng, the inner product {u,v)k can be written as 

(34) 


iu,v)^= Y 

i<jl,---,jd<Ng 


Un,...,33Xju-,3dwY..,j, = 


The Lagrange polynomials Pji,...,jaix) can be used to define d differentiation matrices, defined as 

/ d \ 




/ . 


(35) 


\i=l + l 


Here I is an Ng x Ng identity matrix. can be understood as the discretized differential operator 
di,l < I < d. Similar to Eq. (221, we denote by div a column vector with its entries defined as below 


then div can be expressed in the linear algebra notation as 

div = d\^^v. 

Therefore the inner product (Vu, Vr;)^ can be computed as 


(Vm, Vu)« = I 


(36) 


(37) 




The inner product (u, can be evaluated similar to Eq. (24) as 

d 


(u,u)*,« = + let'll T(u;H)rj 


(38) 


J=i 


with |k| = h‘^. 

In order to evaluate the inner product on the boundary dn, we define d weight vectors corresponding 
to the (d— I) dimensional surface for each dimension I {I = 1,..., d), denoted by w\'^ with the expression 




^ ii,- 

■■,3d y 




,, j; = 1 or ji=Ng, 
l<jl< Ng. 


(39) 


Define Wf = diag 


w 


, then the inner product on the boundary can be expressed as 


iu,v)o^ = 


V 


(40) 
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and 


(Vw • Vi; • n^)a, = 


(41) 


^^=l 


Now we are ready to use the finite dimensional representation of the inner products to evaluate the 
constants d„, a„, b^. 

5.2. Estimation of 

Recall that 


(d^)^ = sup 

«NeVjv(K) 


||Vi;Ar-n„|||^ (VvN-n^,VvM-n^)Q^ 

-= g^p .- 

«iveViv(K) 




{vn, Vn)*,h 


Using Eq. (41) and Eq. (38), we have 


(d^)^ = sup 


vJjMsvn 


Here 


jGVn(k) 'v'n^'>^N 

Ms = 


1=1 


K = + lei'll (lei'll )^. 




(42) 


(43) 


(44) 


Let {^pi{x ),..., ^pn{x)} be a set of basis functions of the finite dimensional space Yn{h)- We denote 
hy ipi{i = 1,..., N) the column vector corresponding to the values of Lpi{x) evaluated at the LGL grid 
points, and denote by 

$ = [(pi,...,(pAr], (45) 

the collection of all column vectors which is an Ng x N matrix. Then for any vector v{x) G the 

corresponding column vector v can be represented as 

V = $c. 


where c is a coefficient vector of size N. Then Eq. (42) can be rewritten as 


(d«,)^ = sup 


c^{^^Ms<S>)c 


cgR« cT{^'^Ks<^)c' 


Eq. (46) can be solved as an eigenvalue problem, 

^^Ms<^c = X^^Ks<^>c, 


(46) 


(47) 


and (d^)^ is equal to the largest eigenvalue A. Since the size of the matrix is N x N and N is 

relatively small, Eq. (47) can be solved as a generalized eigenvalue problem using dense linear algebra. 


5.3. Estimation of 

Recall that 


= sup 


vGH\k), IP 
vJ-Yn (i^) 


{v,v)k 

\*,K, V^H^{k), (^)^)*,K 

vl-Y ]^{k) 


= sup 


then using Eq. (34), Eq. (38) and the density arguments shown above, it can be shown that 

W^Mai; Afo-foo , 

sup 


v^Kv 

ii_LVjv(k) 


(48) 
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where Ma = W^'^\ and K is given in Eq. (44). We can express the orthogonality condition v _L 
in terms of a projection operator Q = 1 — 11^ so that for any v £ Qv _L Yn{k), where I is the 

identity operator. Denoting by 4) as in Eq. (45) the collection of spanning vectors of the space VAr(it), 


then using the Lagrange polynomials corresponding to the LGL grid points as a basis, the projection 
operator 11^ can be expressed as an Ng x Ng matrix 




(49) 


where dt = Therefore the 11^ is a low rank matrix with rank N. The projection 

operator Q and its adjoint operator expressed in the basis of Lagrange polynomials become 


(3 = 7 -$ 4 '^, ( 3 ^ = 7 - 4 '$’ 


Using Eq. (50), the computation of can be simplified as 

2 v'^Q'^M^Qv 

v'^Q^KQv • 


(50) 


(51) 


In other words, corresponds to the largest eigenvalue of the generalized eigenvalue problem 


M^Qv = XQ^KQv. 


(52) 


From a computational point of view, there are two major differences between Eq. (47) and ( |52[) . 
First, the dimension of the matrices in Eq. (47) is iV x TV, and the dimension of the matrices in Eq. (52) 


is Ng X Ng. For 3D simulation, if Ng = 30 then Ng = 27000, and the corresponding eigenvalue problem 
is very costly to solve if M^Q and KQ are treated as dense matrices. Second, the matrix $^774) 


in Eq. (47) is a positive definite matrix since K is positive definite, and the problem (47) can be solved 
directly as a dense generalized eigenvalue problem. On the other hand, Q^KQ is a rank deficient 
matrix with the rank of its kernel being TV. Therefore it can potentially cause a large numerical error 
if Eq. (52) is solved directly as a dense generalized eigenvalue problem. 


In order to overcome the two difficulties mentioned above, we note that for any vector v, the compu¬ 
tational cost for the matrix vector multiplication Qv, Q"^v, M^v, Kv is only proportional to Ng thanks 


to the low rank representation of the operators. Therefore Eq. (52) can be solved using iterative meth¬ 


ods. Another advantage of using iterative methods is that since we only need the largest eigenvalue 
corresponding to Eq. (52), at the fc-th step of the CG iteration we only need to keep three vectors: 


the current approximation of eigenvector the conjugate direction and the residual Even 
though the matrix Q^KQ is singular, the projection onto the 3 dimensional subspace 
is usually well conditioned. In practice we use the Locally Optimal Block Preconditioned Conjugate 
Gradient (LOBPCG) method (with block size equal to 1) for evaluating the largest eigenvalue 
for Eq. (52). It should be noted that since there is no apparent preconditioner that can be applied 


efficiently to solve Eq. (52), the convergence of the largest eigenvalue may be slow. However, we should 
keep in mind that the estimation of aK,bK is only used in the a posteriori error estimator, and only 
low accuracy is needed. In fact is already very accurate in the sense of the preconstant in the 

estimator even if the relative error is 10%. Therefore the slow convergence of the conjugate gradient 
method is compensated by the low accuracy required in the computation of the constants. 

The constant b,.; can be estimated similarly to a^. Recall that 


1,2 

b.., = 


sup 

vGH^{k.), 
vJ-Yn (k) 


#= sup 


{v,v) 


Ok 


V 


{v,v) 


and using the same projection operator Q, b„ can be expressed as 

v'^Q’^M^Qv 


sup 


d v^Q^KQv 


(53) 
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with Mb = Similar to Eq. (521, can be solved as the largest eigenvalue of 

Q^MbQv = XQ^’^KQv. 


(54) 


Eq. (54) can be solved using the same iterative strategy as for obtaining a^. 


6 . Numerical results 

In this section we test the effectiveness of the a posteriori error estimators. The test program is 
written in MATLAB, and all results are obtained on a 2.7 GHz Intel processor with 16 GB memory. 
All numerical results are performed using the symmetric bilinear form {9 = 1). The effectiveness of 
the upper bound and lower bound on the global domain will be justified by comparing |||u — rtjvll 
rj, and by comparing |||m — un\1 and respectively. It should be noted that although our theory does 
not directly predict the effectiveness of the estimator on each local element k, we can measure the local 
effectiveness of the upper and lower bound on each local element k by defining 


C^{k) = 


+ Vf,k + V3,k 

|||w - MmIIU 


C^{k) = 


\u — Un\\ 


(55) 


where the broken energy norm |||u — matUk is defined according to Eq. (|^. 

The numerical results are organized as follows. In section 6.1 we apply the general approach devel¬ 
oped in section to compute the constants for polynomial basis functions, and verify that 

the scaling properties of the numerically computed constants match the analytic results known in the 
literature p4]. In section 6.2 we illustrate the behavior of the upper bound and the lower bound error 


estimates for second order PDEs associated with positive definite operators. We then demonstrate the 
results for indefinite operators in section 6.3 In the a posteriori error estimates of both the upper 


bound and the lower bound, we make the assumption that the non-computable number d“(uAr) can be 
approximated by d^ without significant loss of effectiveness. We justify such treatment in section |6.4| 
by directly calculating d“(itjv) using the numerically computed reference solution. 

Our test problems include both one dimensional (ID) and two dimensional (2D) domains with pe¬ 
riodic boundary conditions. Our non-polynomial basis functions are generated from the adaptive local 
basis (ALB) set 20 in the DG framework. The ALB set was originally proposed to systematically 


reduce the number of basis functions used to solve Kohn-Sham density functional theory calculations, 
and in this section we demonstrate its usage to solve second order linear PDEs. We denote by N the 
number of ALBs per element. For operators in the form of A = —A -|- V with periodic boundary condi¬ 
tion, the basic idea of the ALB set is to use eigenfunctions computed local domains as basis functions 
corresponding to the lowest few eigenvalues. The eigenfunctions are associated with the same operator 
A, but with modified boundary conditions on the local domain. More specifically, in a d-dimensional 
space, for each element k, we form an extended element k consisting of k and its 3“^ — I neighboring 
elements in the sense of periodic boundary condition. On k we solve the eigenvalue problem 

- Atpi + Vifi = Xitpi. (56) 

with periodic boundary condition on dn. This eigenvalue problem can be solved using standard basis set 


such as finite difference, finite elements, or planewaves. Here we solve the local eigenvalue problem (56) 


using planewaves which naturally satisfy periodic boundary conditions. Since this eigenvalue problem 
is solved on a extended element k the computational cost is not large. The collection of eigenfunctions 
(corresponding to lowest N eigenvalues) are restricted from H to k, i.e. 


ipi{x) 


[^i] |„(a;), a; e k; 

0 , otherwise. 


After orthonormalizing locally on each element n and removing the linearly dependent functions, 
the resulting set of orthonormal functions are called the ALB functions. 

Since periodic boundary condition is used on the global domain O, in all the calculations, the reference 
solution, which can be treated as a numerically exact solution, is solved using a planewave basis set with 
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a sufficiently large number of planewaves. The ALB basis set is also computed using a sufficiently large 
number of planewaves on the extended element k. Then a Fourier interpolation procedure is carried out 
from K to the local element k on a Legendre-Gauss-Lobatto (LGL) for accurate numerical integration. 


6.1. Estimating the constants for polynomial basis functions 

Although the main purpose of this paper is to design a posteriori error estimator for non-polynomial 
basis functions, the computational strategies discussed in section can be applied to polynomial func¬ 
tions as well. Let k = [0, and Yn{p]k) = span{S — P} be the space spanned 

by polynomials with degree less than or equal to p. Then the asymptotic scaling of a,.j,bK,dK with 
respect to h and p is known 


a 


2 

K 


— h2 


h 

7 

P 



(57) 


These results are asymptotically correct as p —>■ oo, and we will show that the strategy discussed in 
section leads to the same asymptotic result, but the result is more accurate in the pre-asymptotic 
regime due to the explicit computation of the constants. 

From numerical point of view, the scaling with respect to h is naturally satisfied. To verify this, 
we can simply consider a reference element K\h=i = [0,1]'^ and scale the weight matrix and the 
differentiation matrix d\‘^^ accordingly. The technique is the same as that used in 


24 


We now directly verify the scaling with respect to p in Fig. using the algorithms presented in 
section The LGL grid sizes for ID, 2D and 3D calculation are chosen to be 100, 100 x 100, and 
50 X 50 X 50, respectively. The largest degree of polynomials is 64 for ID and 2D, and is 16 for the 
3D case. Note that in the 3D case, the dimension of Vat( p = 16; k) is already 969. Fig. (a) shows 
the behavior of a^, which asymptotically agrees with the 1/p^ scaling. It is interesting to see that the 
computed can be approximated by where the constant C is around 0.1. The recovery of the 
constant indicates that the numerically computed constant a^ can offer a sharper estimator even for 
the standard hp-refinement. Similarly Fig. [^(b) shows that asymptotically scales as 1/p for 2D and 
3D simulation. The ID case is not shown in the picture, since the numerical value of b^ is already 
as small as 10“^° for p = 2. This can be interpreted from Proposition 7.1 in the appendix. Finally, 
direct computation in Fig. (c) shows that asymptotically scales as p^ for all dimensions. Again, 
the computed constant differs from the asymptotic scaling in the pre-asymptotic regime, indicating 
that the numerically computed constant should be sharper for low order polynomials (p < 4). 
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p-refinement 


p-refinement 

10 “ -^ 1 ^- 



p-refinement 



(c) 


Figure 1. Numerically computed constants a^,b^,d^ with respect to the polynomial 
degrees p in ID, 2D and 3D. 


6.2. Positive definite operators 

We first demonstrate the effectiveness of the a posteriori error estimates for a positive definite operator 
on a ID domain = (0, 27r), using the ALB set as non-polynomial basis functions. Due to the periodic 
boundary condition, we choose V{x) = 0.01 so that the operator A = —A -|- F is non-singular and 
positive definite. The right hand side is chosen to be f{x) = sin( 6 a;) which is periodic on P. In the 
ALB computation, the domain is partitioned into 7 elements, as indicated by black dashed lines. Fig. 
shows solution u to Eq. (14) and the point-wise error u — using A = 11 ALBs per element. 

Fig. i (a) shows the absolute error in the energy norm, the upper bound and lower bound estimates 
as the number of ALBs per element N increases from 3 to 15. The relative error can be deduced by 
comparing Fig. (a) and Fig. (a). We find that the computed 77 and ^ are indeed upper and lower 
bounds of the true error |||m — un\1 for all N across a wide range of accuracy (from 10“^ to 10“®). 
It also appears that the lower bound estimator ^ follows the true error more closely than the upper 
bound estimator p. Fig. |^(b) and (c) illustrate the local effectiveness C^(«;) and C^{k) for each element 
K. Though not guaranteed by our theory, we observe that and are upper and lower bounds for 
|||w — mtvIIU fo'' each element k, respectively. The effectiveness as measured by Cjj{k) and C^{k) depends 
only weakly on the number of adaptive local basis functions, or the accuracy of the numerical solution. 

Our next example is to solve a 2D problem with O = (0, 27r) x (0, 27r). Again we choose V{x, y) = 0.01 
so that A = — A-l-F is non-singular and positive definite. The right hand side is f{x, y) = cos(3x) cos{y), 
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(a) (b) 

Figure 2. (a) The reference solution u{x) corresponding to V{x) = 0.01 and the right 
hand side f{x) = sin(6a:). (b) Point-wise error between the reference solution u{x) and 
the numerical solution un{x) calculated using the ALB set with 7 elements and iV = 11 
basis functions per element. The domain is partitioned into 7 elements indicated by 
black dashed lines. 


which satishes the periodic boundary condition. Fig. shows the reference solution u to Eq. (141 and 
the point-wise error u — un using = 31 ALBs per element. In the ALB computation, the domain is 
partitioned into 5x5 elements, indicated by black dashed lines. 

Fig. El (a) shows the error in the energy norm, the computed upper bound and the lower bound as 
the number of ALBs per element N increases from 11 to 41. Both the computed upper and the lower 
bound estimates are effective for all calculations. Fig. (b)-(d) illustrates the local effectiveness of the 
upper and lower bound estimates for the two extreme cases A = 11 and A = 41, and the estimator 
77 k and are effective for all elements, and the effectiveness depends weakly on the number of basis 
functions per element. 
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(a) (b) 



Figure 3. (a) Global error and the upper/lower bound estimator for V{x) = 0.01 and 
f{x) = sin(6a;). (b) Local effectiveness of the upper bound characterized by C^{k) for 
each element, (c) Local effectiveness of the lower bound characterized by C^{k) for 
each element. 
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Figure 4. (a) The reference solution u{x,y) corresponding to V{x,y) = 0.01 and 
f{x,y) = cos(3a;) cos(i/). (b) Point-wise error between the reference solution u{x,y) 
and the numerical solution uiq{x,y) calculated using the ALB set with 5x5 elements 
and = 31 basis functions per element. 
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Figure 5. (a) Global error and the upper/lower bound estimator for V{x,y) = 0.01 
and f{x,y) = cos(3a;) cos(?/). (b) Local effectiveness of the upper bound characterized 
by Cri in each element for = 11. (c) Local effectiveness of the upper bound charac¬ 
terized by Crj in each element for TV = 41. (d) Local effectiveness of the lower bound 
characterized by in each element for N = 11. (e) Local effectiveness of the lower 
bound characterized by in each element for N = 41. 
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6.3. Indefinite operators 


We now demonstrate the effectiveness of the upper and lower bound estimates for indefinite operators. 
We start from a ID example on a domain = (0,27r) with periodic boundary conditions. The potential 
function V(x) is given by the sum of three Gaussians with negative magnitude, as shown in Fig. 
(a). The operator A = —A + V has 3 negative eigenvalues and is indefinite. The right hand side is 
f{x) = sin(6a;). The domain is partitioned into 7 elements for the ALB calculation. Fig. (b) shows 
the reference solution u to Eq. (141, and Fig. (c) shows the point-wise error u — un using = 11 
ALBs per element. 




(a) (b) 



Figure 6 . (a) The potential V{x) given by the sum of three Gaussians with negative 
magnitude, (b) The reference solution u{x) corresponding to the potential V(x) in (a) 
and the right hand side f{x) = sin(6a;). (c) Point-wise error between the reference 
solution u{x) and the numerical solution mat (a;) calculated using the ALB set with 7 
elements and A^ = 11 basis functions per element. 

Fig. 0 (a) shows the error in the energy norm, the computed upper and lower bound estimates as 
the number of ALBs per element N increases from 3 to 15. Similar to Fig. the computed rj and 
^ are upper and lower bounds for the true error |||m — unW for all N across a wide range of accuracy. 
Furthermore, the computed ^ is always a lower bound of |||u — utvI from N = 3 to N = 15. This is 
guaranteed by the property of the lower bound in Proposition |4.5| 

We should note that when the number of basis functions is very small {N = 3), the accuracy is low 
and the ALB approximation is in its pre-asymptotic regime. In such case, the upper bound is very close 
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to the true error. In fact as indicated by Theorem 4.1 rj may not even be a rigorous upper bound for 
highly indefinite operators with very few basis functions. 




Figure 7. (a) Global error and the upper/lower bound estimator for V{x) given in 
Fig.§ (a) and f{x) = sin(6x). (b) Local effectiveness of the upper bound characterized 
by Cn in each element, (c) Local effectiveness of the lower bound characterized by Crj 
in each element. 


Our final examples are two indefinite problems on a 2D domain O = (0,27r) x (0, 27r). The first 
problem is a homogeneous Helmholtz equation with V{x,y) = —16.5 and the operator A = —A + V 
has 49 negative eigenvalues. The right hand side is 

f{x, y) = exp(-2(a; - iif - 2{y - irf), (58) 


which is a Gaussian located at the center of O. The second problem is that V is given by the sum of 
four Gaussians with negative magnitude, as illustrated in Fig. 10 (a). The operator A = — A + F has 26 
negative eigenvalues. The right hand side is chosen to be f{x, y) = cos(3x) cos{y) satisfying the periodic 
boundary condition. For the first problem, Fig. (b) shows the reference solution u to Eq. (14) and 
Fig. § (c) shows the point-wise error u — un using iV = 31 ALBs per element. In the ALB computation, 
the domain is partitioned into 5x5 elements, indicated by black dashed lines. Similarly for the second 
problem. Fig. 10 shows solution u to Eq. (14) and the point-wise error u — un using iV = 31 ALBs per 
element. 

Fig.|^ (a)-(e) illustrates the global and local effectiveness of the upper and lower bound estimates for 
the Helmholtz problem, as the number of ALBs per element N increases from 21 to 51. Compared to 
the positive definite case in Fig.[^ the true error is larger using a comparable number of basis functions, 
reflecting that the Helmholtz equation is more difficult to solve. Nonetheless, rj and ^ provide effective 
bounds for the true error in all cases. Similar results can be found for the indefinite example with 
negative Gaussian potentials in Fig. 11 (a)-(e). In all calculations, the computed lower bound estimator 
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remains a lower bound for the true error. In particular, the estimators still hold quite tightly in the 
pre-asymptotic regime {N = 11) where the ALB approximation is crude and has large numerical error. 
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Figure 8. (a) The reference solution u{x,y) corresponding to V{x,y) = —16.5 and 


/(x, y) in Eq. (58), which is a Gaussian localized at the center of fl. (b) Point-wise error 


between the reference solution u{x,y) and the numerical solution UN{x,y) calculated 
using the ALB set with 5x5 elements and A^ = 31 basis functions per element. 
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Figure 9. (a) Global error and the upper/lower bound estimator for V{x,y) = —16.5 
and f{x,y) in Eq. (58), which is a Gaussian localized at the center of 12. (b) Local 
effectiveness of the upper bound characterized by in each element for N = 21. 
(c) Local effectiveness of the upper bound characterized by in each element for 
= 51. (d) Local effectiveness of the lower bound characterized by in each element 
for N = 21. (e) Local effectiveness of the lower bound characterized by in each 
element for = 51. 
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Figure 10. (a) The potential V{x,y) four Gaussians with negative magnitude, (b) 
Solution u(x,y) corresponding to V{x,y) given in (a) and f{x,y) = cos(3a;) cos(?/). 
(c) Point-wise error between the reference solution u{x, y) and the numerical solution 
un{x, y) calculated using the ALB set with 5x5 elements and = 31 basis functions 
per element. 
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Figure 11. (a) Global error and the upper/lower bound estimator for V{x,y) given 
in Fig. [^(a) and f{x,y) = cos(3ai) cos( 2 /). (b) Local effectiveness of the upper bound 
characterized by Cjj in each element for N = 11. (c) Local effectiveness of the upper 
bound characterized by in each element for N = 41. (d) Local effectiveness of the 
lower bound characterized by in each element for = 11. (e) Local effectiveness of 
the lower bound characterized by in each element for = 41. 
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6.4. Justification of the treatment of d“(nAr) 

In the numerical computation of the upper and lower bound estimates, we approximated the non- 
computable constant d^(ujY) by the computable constant d^. Below we provide numerical justification 
of such approximation by direct computation of d“(itjv) via the reference solution. We compare with 
d^ and since these three terms appear together in in Eq- (121. 

Fig. 12 (a) and (b) compare d“(u 7 v), d,^ and for the positive definite and the indefinite ID 
examples, respectively. We observe that the magnitude of d“(uAr) is comparable to that of d,.;. is 
much smaller compared to d“(rtjv) and d^. This is a direct consequence of Proposition 7.1 which states 
that b„ is in general very small for ID systems. 



Figure 12. Comparison of dJ^(ujY), d,^ and b^yK for (a) the positive definite case with 
V{x) = 0.01 with N = 7. (b) the indehnite case with V{x) given in Fig. (a) with 
N = 7. 


Fig. 13 compare d“(ujv), d^ and b„y„ for the positive definite case V = 0.01, the indefinite case 
V = —16.5, and the indefinite case with V given by the sum of negative Gaussians in Fig. (a). In 
all cases, the magnitude of d“(uAr) is comparable to that of d^. Furthermore, both d“(ujv) and d^ are 
much smaller compared to PkIk- Therefore the effectiveness of the estimator remains unchanged even 
if d“(u 7 v) is neglected. We expect similar results can be observed for systems of higher dimensionality. 

Finally we provide a second justification by comparing the total contribution of the jump term in 
the upper bound estimator 




= E 




and the total contribution of the jump term in the energy norm 


Ej = Y.Y\\M\ 


2 

Ok' 


This is given in TableIt shows that the approximation ~ d^ does not lead to underestimation 

of the jump term, which is consistent with the observation in Fig. |12|and 13 
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(a)y = 0.01 


(h)V = 0.01 


Y = 21 



(c)V = 0.01 





(g)V' Gaussian 


(h)V' Gaussian 


{i)V Gaussian 


Figure 13. Comparison of d“(ujv), and for 2D test problems for (a-c) the 
positive definite case V = 0.01 (d-f) the indefinite case V = —16.5 (g-i) the indefinite 
case with V given by the sum of negative Gaussians in Fig. 10 (a). 


Problem 

N 

Ej 

2 

Vj 

ID F = 0.01 

7 

2.0179 X 10-« 

2.0182 X 10-3 

2D F = 0.01 

21 

1.2030 X lO"'^ 

9.1593 X 10-3 

ID Gaussian 

11 

6.4687 X 10"“ 

6.4697 X 10-“ 

2D F = -16.5 

31 

4.7352 X 10-3 

5.6649 X 10-2 

2D Gaussian 

21 

1.6226 X 10-3 

2.8348 X 10-2 


Table 1. Comparison of the total contribution of the jump term in the estimator rij, 
and the total contribution of the jump term in the energy error Ej. 


















































































32 


TITLE WILL BE SET BY THE PUBLISHER 


7. Conclusion and future work 

We present the first systematic work for deriving a posteriori error estimates for general non¬ 
polynomial basis functions in a interior penalty discontinuous Galerkin (DG) formulation for solving 
second order linear PDEs. The estimates not only serve to quantify the error sharply for a given com¬ 
putation, but also can lead an adaptive algorithm to refine the elements non-uniformly by adding (or 
even removing/coarsening) basis functions to certain elements. This allows a best approximation for 
a given number of degrees of freedom in order to reduce the computing time even when relatively few 
degrees of freedom are employed. A non-uniform distribution of the number of local basis functions is 
in this case mandatory to develop powerful solvers, in particular when inhomogeneous data of the PDE 
is involved. It turns out that the standard polynomial hp DG-method may benefit from this analysis 
as it involves numerically computed constants. 

Our analysis requires the exact solution to lie in H^{k) for each element k which may seem limiting 
when dealing with a posteriori estimates for Poisson’s equation as a uniform refinement leads to optimal 
convergence rates in the asymptotic limit. We remark that despite the above asymptotic reasoning there 
are numerous cases where the a posteriori analysis for regular functions is still interesting, for example if 
the PDE involves a strong small-scale character (but still being smooth) either due to strongly oscillating 
material coefficients or a wave-like character of the underlying PDE (Helmholtz equation for instance). 
Or, if the data of the PDE and thus the solution as well has an inhomogeneous character so that a 
uniform refinement involves too many degrees of freedom. In this case, combining the estimates with 
an adaptive algorithm as outlined above will result in an optimal balance of degrees of freedom per 
element. 

Our framework for developing explicitly computable constants for a posteriori error estimates are 
not limited to second order PDEs, nor it is necessarily limited to discontinuous Galerkin framework. In 
a forthcoming publication we will demonstrate the method for eigenvalue problems. It is also possible 
to generalize the method to multiscale methods and reduced basis methods. 
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Appendix 

Proposition 7.1. Let k = [a, b] be a ID element and VAr(_p; k) = span{x^ , j < p} be the function space 
spanned by polynomials with degree less than or equal to p. Then Vp > 2, = 0. 

Proof. Define c = (a -I- b)/2. For any v € H^{k), v T 'Vn{p] k) with p > 2, we have 
(i;,l)*,K=0, (v, (a; - c))*,„ = 0, (i;, (x - c)^)*,„ = 0. 

Using the dehnition of the inner product (•, •)*_„ 

pb pb pb 

/ v{x)dx = 0, / v'{x)dx = 0, / v'{x){x — c)dx = 0. 

J a J a J a 

With integration by parts, we have v{a) = v{b) = 0. Therefore ||r|| 3 K = 0. Using the definition of b„ 
we obtain b^ = 0. □ 
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